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ABSTRACT 

We investigate the vorticity of the IGM velocity field on large scales with 
cosmological hydrodynamic simulation of the concordance model of ACDM. We 
show that the vorticity field is significantly increasing with time as it can effec- 
tively be generated by shocks and complex structures in the IGM. Therefore, the 
vorticity field is an effective tool to reveal the nonlinear behavior of the IGM, 
especially the formation and evolution of turbulence in the IGM. We find that 
the vorticity field does not follow the filaments and sheets structures of under- 
lying dark matter density field and shows highly non- Gaussian and intermittent 
features. The power spectrum of the vorticity field is then used to measure the 
development of turbulence in Fourier space. We show that the relation between 
the power spectra of vorticity and velocity fields is perfectly in agreement with the 
prediction of a fully developed homogeneous and isotropic turbulence in the scale 
range from 0.2 to about 3h~ l Mpc at z ~ 0. This indicates that cosmic baryonic 
field is in the state of fully developed turbulence on scales less than about 3 hr x 
Mpc. The random field of the turbulent fluid yields turbulent pressure to prevent 
the gravitational collapsing of the IGM. The vorticity and turbulent pressure are 
strong inside and even outside of high density regions. In IGM regions with 10 
times mean overdensity, the turbulent pressure can be on an average equivalent 
to the thermal pressure of the baryonic gas with a temperature of 1.0 x 10 5 
K. Thus, the fully developed turbulence would prevent the baryons in the IGM 
from falling into the gravitational well of dark matter halos. Moreover, turbu- 
lent pressure essentially is dynamical and non-thermal, which makes it different 
from pre-heating mechanism as it does not affect the thermal state and ionizing 
process of hydrogen in the IGM. 

Subject headings: cosmology: theory - intergalactic medium - large-scale struc- 
ture of the universe - methods: numerical 
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1. Introduction 

Gravity is curl-free in nature, therefore it is unable to trigger vorticity within the velocity 
field of a cosmic flow. On the linear order of cosmological perturbation theory, the vorticity 
will inevitably decay due to the expansion of the universe. The linear velocity fields of the 
cosmic flow should be irrotational. In the nonlinear regime of clustering, vorticity can be 
generated in the collisionless dark matter field when multi-streaming occurs at shell crossing 
(Binney, 1974, Pichon & Bernardeau, 1999). However, there is no way to directly map 
the vorticity of the dark matter field to the baryon field, most of which is the intergalactic 
medium (IGM). 

In the context of fluid dynamics, vorticity can be generated if the gradient of the mass 
density and the pressure gradient of cosmic flow are not aligned (Landau & Lifshitz 1987). 
Namely, vorticity results from the complex structures of fluid flow like curved shocks. Re- 
cently, it has been revealed that in the nonlinear regime, the cosmic baryon fluid at low 
redshift does contain such complex structures (e.g. He et al 2004). Therefore, one expects 
that vorticity would be present and evolve extensively in the cosmic baryonic field. The 
vorticity of the intracluster medium (ICM) has been studied in topics related to possible 
mechanism of generating magnetic field of galaxies or clusters (Davis & Widrow 2000; Ryu, 
et al 2008). Although these works show that the vorticity can form in the ICM, the formation 
and evolution of vorticity in the IGM is still unknown. 

In addition, no studies have been done on the relation between vorticity and turbulence 
in a cosmic baryon fluid. Actually, vortices generally are considered a fundamental ingredient 
of turbulence and the fluctuations of the vorticity field is an important indicator to describe 
the turbulence of fluid (e.g. Batchelor, 1959, Schmidt 2007). On the other hand, the study 
of the turbulence of cosmic fluid on large scales has seen a lot of progress in recent years. 
The fluctuations of the velocity field of the baryon fluid beyond the Jeans length is shown to 
be extremely well described by the She-Leveque (SL) scaling (He, et al 2006), which is the 
generalized scaling of the classical Kolmogorov's 5/3-law of fully developed turbulence (She 
& Leveque 1994). The non-Gaussian features of the density field in baryon flows are found to 
be in good agreement with the log-Poisson cascade (Liu & Fang 2008), which characterizes 
statistically the hierarchical structure in fully developed turbulence (Dubrulle 1994; She & 
Waymire 1995; Benzi et al. 1996). Observationally, the intermittence of Lya transmitted 
flux of QSO absorption spectrum can also be well explained in terms of log-Poisson hierarchy 
cascade (Lu et al. 2009). These results suggest that the dynamical behavior of the IGM is 
similar to a fully developed turbulence in inertial ranges. Therefore, it would be worthwhile 
to investigate the vorticity fields of the turbulent cosmic fluid. 

An important problem related to the vorticity fields of a turbulent fluid is the impact 
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of the turbulent pressure on the clustering of cosmic fluid. It is well known that the random 
velocity field of a turbulent fluid will play a similar role as thermal pressure and prevent 
the gravitational collapse in such a fluid (Chandrasekhar, 1951; Bonazzola et al. 1992). In 
the ICM, this effect has been studied with hydrodynamic simulations (Dolag et al 2005; 
Iapichino & Niemeyer 2008; Cassano 2009, and reference therein). However, in these works 
the turbulent pressure is directly identified with the RMS baryon velocity. This identification 
may be reasonable for the ICM; however, it would be a poor relation on scales larger than 
clusters, as the RMS baryon velocity cannot separate the velocity fluctuations due to bulk 
motion from that of turbulence. Obviously, the bulk motion is not going to prevent gravita- 
tional collapsing. Since the dynamical equation of vorticity is free from gravity, the vorticity 
field provides an effective method to pick up the velocity fluctuations within a turbulent 
flow. The power spectrum of the vorticity field yields a measurement on the scale of velocity 
fluctuations where turbulence is fully developed. Using this method, we can estimate the 
turbulent pressure in the IGM and hence study its effect on gravitational clustering. 

We will investigate these problems with cosmological hydrodynamic simulation samples 
of the concordance model of ACDM. In §2, we present the equations governing the dynamics 
of vorticity and rate of strain field. §3 gives a brief description of the cosmological hydrody- 
namic simulation of the ACDM model. In §4 we discuss the statistical properties of vorticity 
on large scales. The nonthermal pressure of turbulent fluid and its effects on clustering of the 
IGM are addressed in §5. We summarize the basic results of the paper and give concluding 
remarks in §6. Mathematical equations are given in the Appendix. 



The dynamics of a fluid is conventionally governed by a set of equations for velocity 
and density fields Vi(t,r), p(t,r) (Appendix §A.l). An alternative way is to replace the 
velocity field by their spatial derivatives diVj. The velocity derivative tensor d{Vj can be 
decomposed into a symmetric component SV,- = (l/2)(<9jfj + djVi) and an antisymmetric 
component (l/2)(9ji;_ 7 - — djVi) (Landau & Lifshitz 1987). The former is the rate of strain and 
the latter is the vorticity vector Ui = eijkdjVk, or uj = V x v, where is the Levi-Civita 
antisymmetric symbol. For a cosmic baryon fluid (IGM), the dynamical equation of vorticity 
uj can be derived from the Euler equation as (Appendix §A.l and §A.2) 



2. Theoretical Background 



2.1. Dynamical Equation of Vorticity 



Du 




(1) 



Dt 
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where p is the pressure of the IGM, a(t) is the cosmic factor, d = diVi is the divergence of 
the velocity field, and the vector [S • u]i = SyOUj. Defining a scalar field as u> = |cJ|, the 
dynamical equation of oj is then 



Duj o 1 „ 1 
— — = d t u; + -v • Vw = - 

Dt a a 



where £ — uj/uu, and a = £ • (£ • V)v. 



aa> — da; H — (Vp x Vp) — acu 



(2) 



An essential feature of both eqs.(l) and (2) is that they are free from the gravity of 
mass fields, therefore, the gravitational field of both dark matter and the IGM cannot be a 
source of the vorticity. Obviously, in the linear regime, only the last term of eqs.(l) and (2) 
survives. This term is from the cosmic expansion, and makes the vorticity decaying as a -1 . 
Thus, the vorticity of the IGM is reasonably negligible in the linear regime. 

Equations (1) and (2) show that if the initial vorticity is zero, the vorticity will stay at 
zero in the nonlinear regime, provided that the term (1/p 2 ) Vp x Vp is zero. This term, called 
baroclinity, characterizes the degree to which the gradient of pressure, Vp, is not parallel to 
the gradient of density, Vp. 

If the pressure of a baryon gas is a single-variable function of density, e.g. there exists 
a determined relation for the equation of state p = p(p), the vector Vp would be parallel 
to Vp, and then (l/p 2 )Vp x Vp = 0. Therefore, vorticity cannot be generated even in 
the nonlinear regime until the single-variable function or determined relation for p = p(p) 
is violated. Physically, once multi-streaming and turbulent flows have developed, complex 
structures, like curved shocks, will lead to a deviation of the direction of Vp from that of 
Vp. In this case, the p — p relation cannot be simply given by an single-variable function 
equation as p = p(p) (He et al 2004) and the baroclinity will no longer be zero. 

The term S • uj on the right hand side of eq.(l) accounts for stretch of vortices drived 
by strain. The vorticity will be either amplified or attenuated by this term. Actually, this 
point can be easily seen with eq.(2). If the coefficient a is larger than zero, i.e. £ is in the 
direction of the eigenvector of tensor djVi with positive eigenvalue, the vorticity will grow at 
the rate of au. Otherwise it would be attenuated. The term — du stands for expansion or 
contraction of vortices caused by the compressibility of baryon. Since divergence d = djVj 
is generally negative in regions of clustering, the term — du will lead to an amplification of 
vorticity in overdense regions. 
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2.2. Vorticity Effect on IGM Clustering 

The effect of vorticity on the IGM clustering can be seen from the dynamic equation of 
divergence d, which is an indicator of clustering. The equation reads (Appendix §A.3) 

E* = fifcd+IyVd (3) 

Dt a 



\u 2 - S^S^ - -V 2 p + 4(Vp) • (Vp) - —(ptot - Po) - ad 



where p to t is the total mass density including both CDM and baryon and po is its mean 
value. A negative d corresponds to a convergent flow (clustering), while a positive d means 
a divergent flow. As in the equation (2) for vorticity, there is a term — ad coming from the 
cosmic expansion that leads to dilution of d. However, different from eqs.(l) and (2), the 
gravity effect —4nG(p to t — Po)/a, Po ac ts as a source term in the divergence equation. This 
term leads to clustering in regions with p tot > po, and anti-clustering for p tot < po- 

The term (Vp) • (Vp)/p 2 will be nonzero even when the IGM is barotropic, or the 
density-pressure relation is a power law p oc p 7 and 7 > 0. The ratio of this term to the 
gravity is roughly ~ (t infaU /t sound ) 2 , where t infaU ~ (Gp)~ 1/2 , and t sound ~ l/c s with the 
typical scale of density variation I ~ (Vp/p) -1 and the speed of sound c s ~ (Vp/Vp) 1 / 2 
. The value of this ratio defines roughly the Jeans criterion for gravitational instability. 
In addition, the pressure term — V 2 p is compatible with (Vp) • (Vp)/p 2 and is likely to be 
positive in overdense clustering regions. Hence, these two terms are from thermal pressure 
to resist upon gravitational collapse. 

Finally, we examine the effect of the first two terms, the strain S^Sy and the vorticity 
|w 2 , on the right hand side of eq.(3). For simplicity, we consider an incompressible fluid in 
the absence of gravity. In this case, eq.(3) simplifies to 

V 2 p = -p (s^Sy - \^ . (4) 

This is a typical Poisson equation for a scalar field of the pressure p. Taking the similarity 
with the field equations in electrostatics, the term on the right hand side of eq.(4), = 
p[SijSij — l/2u 2 ], mimics the "charge" of a pressure field. A positive "charge" produces an 
attraction force that tends to drive overdense charge halos while a negative "charge" yields 
a repulsive force that smear out the charge accumulation. Back to the IGM flow, £} plays 
the role of nonthermal pressure of turbulence (Chandrasekhar 1951, a, b; Bonnazzola et al 
1987). In regions with £3 < 0, the turbulent pressure will prevent the IGM clustering. The 
sign of H is actually determined by levels to which the turbulence has developed(§5.2). 
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3. Numerical Method 

To model the flow patterns of the IGM and dark matter fields, we use the WIGEON 
code, which is a cosmological hydrodynamic/N-body code based on the fifth-order WENO 
(weighted essentially non-oscillatory) scheme (Fenget al. 2004). The WENO scheme uses the 
idea of adaptive stencils in the reconstruction procedure based on the local smoothness of the 
numerical solution to automatically achieve high order accuracy and non-oscillatory property 
near discontinuities. Specifically, WENO adopts a convex combination of all the candidate 
stencils, each being assigned a nonlinear weight which depends on the local smoothness of 
the numerical solution based on that stencil (Shu, 1998, 1999). For more details, one can 
refer to Appendix A. 4. 

The WENO scheme has been successfully applied to hydrodynamic problems containing 
turbulence (Zhang et al 2008), shocks and complex structures, such as shock- vortex interac- 
tion (Zhang et al 2009), interacting blast waves (Liang & Chen 1999; Balsara & Shu 2000), 
Rayleigh- Taylor instability (Shi, Zhang & Shu 2003). The WENO scheme has also been used 
to simulate astrophysical flows, including stellar atmospheres (del Zanna, Velli & Londrillo 
1998), high Reynolds number compressible flows with supernova (Zhang et al. 2003), and 
high Mach number astrophysical jets (Carrillo et al. 2003). In the context of cosmological 
applications, the WENO scheme has been proved to be especially adept at handling the 
Burgers' equation, a simplification of Navier-Stokes equation,typically for modeling shocks 
and turbulent flows (Shu 1999). This code has also been successfully applied to reveal the 
turbulence behavior of the IGM (He et al 2006, Liu & Fang 2008, Lu et al 2009). 

We evolve the simulation in the concordance model of a LCDM universe specified by the 
cosmological parameters (fi m , Sl A , h, a 8 , Sl b , n s , z re ) = (0.274, 0.726, 0.705, 0.812, 0.0456, 0.96, 11.0) 
(Komatsu et al., 2009). The simulation is performed in a periodic cubic box of size of 25 h^ 1 
Mpc with a 512 3 grid and an equal number of dark matter particles, which have mass reso- 
lutions 1.04 x 1O 7 M . To test the convergence, we also run a low-resolution simulation with 
a 256 3 grid and an equal number of dark matter particles in the same box. Radiative cooling 
and heating are modeled using the primordial composition (X = 0.76, y = 0.24) and calcu- 
lated as in Theuns et al.(1998). A uniform UV background of ionizing photons is switched 
on at z re . Processes such as star formation and feedback due to stars, galaxies and active 
galactic nuclei(AGN) are not included in our simulation. The simulations start at redshift 
z = 99, and the snapshots are outputted at redshifts z = 11.0, 6.0, 4.0, 3.0, 2.0, 1.0, 0.5, 0.0. 

The tensor diVj of samples is then calculated by using a four-point finite-difference 
approximation at the same grid that is used in the simulation. For example, the partial 
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derivatives of d y v x at grid l,m, n is given by 

2 1 

dyV x (l,m,n) = -[v x (l,m + l,n) -v x (l,m- l,n)\ - — [v x (l,m + 2,n) - v x (l,m- 2,n)\. (5) 

Once all the partial derivatives of three velocity components are generated, one can produce 
the fields of vorticity and the rate of strain of these samples. 



4. Basic Properties of the IGM Vorticity 
4.1. Configuration of the Vorticity Fields 
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Fig. 1. — 3-D distribution of density fields of dark matter(left) and baryon(right) in a periodic 
box size of 25/i" 1 Mpc with a 512 3 grid. 

Figure 1 visualizes 3-D density distributions of the dark matter (left) and the baryon 
(right) respectively. Figure 2 gives the 3-D distributions of scalar field ut,t is the cosmic 
time, at redshifts z = 4, 2 and 0. The dimensionless quantity cut is to characterize the 
typical number of eddy turnovers of vorticity within the cosmic time. Figure 2 shows a 
strong evolution of the intensity of vorticity with redshifts. The vorticity has not been well 
developed until redshifts z ~ 4, but becomes significant at z ~ 2 which marks approximately 
the onset of shocks and complex structures developed on the cosmic scales, and matches with 
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the typical formation history of dark matter halos of galactic clusters (e.g. Bahcall & Fan 
1998). 



lg(l»lt) z= 2 



lg(l<alt) z= 




Fig. 2. — 3-D distribution of dimensionless scalar field Lot that refers to the number of 
turnovers of vortices within the cosmic time, where t is the cosmic time, at redshifts z = 4 
(left), 2 (middle) and (right). 

The density fields (Figure 1) display the typical sheets-filaments-knot structures on the 
cosmic scale. However, the spatial configuration of the vorticity field looks quite different: 
it does not show any sheet-like or filamentary structure, instead, looks like clouds with the 
comparable sizes of clusters. Although the vorticity field is not associated with the fine 
structures of the underlying density fields, the cloudy structures are most likely to occur 
around overdense regions on large scales, and thus the vorticity field can be used to pick up 
the coherent structures. These features can be more clearly illustrated by 2-D distributions 
of cut and vector projection of lo shown in Figure 3, where a slide of 25 x 25 x 0.1 h~ 3 Mpc 3 
is taken. The distribution of lo shows a similar spatial pattern as the scalar quantity ut. 

It is recalled that, in a incompressible fluid, the evolution of vorticity is driven domi- 
nantly by the amplification of the strain rate term S-lo [eq.(l)], which tends to stretch the vor- 
tical motion (Tanaka, & Kida, 1993; Constantin et al 1995) and produce filamentary(tubes) 
and/or sheetlike structure in the vorticity field (e.g. She et al 1990). As mentioned in §2.2, 
the strongest stretching is in the direction of £, which is parallel to the eigenvector of tensor 
djVi with large positive eigenvalue. This mechanism distorts the vorticity field and forms a 
tube-like network in the spatial configuration. 

The IGM, however, is compressible as a fluid. The amplification due to the strain rate 
S ■ lo will be largely canceled by the term —dto [eq.(2)], which results in a strong attenua- 
tion of vorticity in the direction parallel to the eigenvector of the tensor djVi with positive 
eigenvalues. Consequently, the vector field lo does not show any clear sheetlike-filamentary 



-9 - 



'■gfrb) z= 



'g(PHlni) z= 
— I — I — I — I — I — 



1.0 

{3,8 
0,6 



4.Q -3.0 -2.0 -1.0 0.0 1.0 , 2.0 3.0 4.0 -4.Q -3.0 -2.0 -1.0 0.0 1.0 , 2.0 3.0 4.0 



7":'. <v 



t / 
-4 



\ .1 



j . -_ 



lh*- ^ : is! 



0,2 



0-0 P*-' i 



'A 



■ I 



1.0 
0,B 

V 

0,6 



0,2 



0,0 kiL 



1 



■7* J •AJ^N" 



y - : l 

- % s £ : 



■it- '; - 



- 



D.Q 0.2 0.4 0.6 O.B 1.Q D.Q 0.2 0.4 0.6 O.B 1.0 



Iq(blt) z= 




lgll/p 2 VpxVp| z = □ 

1 — 1 — i — 1 — i — 1 — i — 1 — 



-6. 0-5.0 -4.0 -3.0 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 S .Q 
1.0 I — ' — : — ' — I ' ,| ""' — ' — I — ' — ' — ' — I — ' — ' 




4 



■ 



i 



0.2 



D.4 0.6 0.8 1.0 



Fig. 3. — Vorticity in a slide of 25 x 25 x 0.1 /i~ 3 Mpc 3 . The top two plots give vector fields 
of cu against background of baryon density(top left) and dark matter density (top right). The 
bottom left plot presents ut in this slide while the baroclinity field is given by the bottom 
right panel. 



- 10 - 



structure. 

Since vorticity is mostly attributed to the baroclinity (l/p 2 )Vp x Vp, the distribution 
of vorticity should be determined by the distribution of baroclinity. In figure 3, we also 
present the baroclinity field in the same slide as that of ut. Clearly, both of them show alike 
structures. It is noted that, similar to the vorticity, the baroclinity can be strong even at 
low density regions, as shocks and complex structures can be formed there (He et al 2004). 

Nevertheless, there does not exist a linear mapping between the ut and the baroclinity. 
This is because the term \au — du\ is sometimes comparable with the baroclinity |(l/p 2 )£ • 
(Vp x Vp)|. Figure 4 gives a cell-by-cell comparison between \au — du\ and |(l/p 2 )£ ■ (Vp x 
Vp)|. In average, the intensity of these two sources are almost of the same order. Thus, the 
amplification and stretching by the rate-of-strain and divergence cannot be ignored. It leads 
to the mapping between the ut and the baroclinity field deviating from a linear one. 




BET 3 10"* 10" 1 10° 10 1 10 s 10 3 
l1//i 3 V/ixVpl 



Fig. 4. — A cell by cell comparison between the term \au — du\, accounting for stretch in 
addition to expansion or contraction of vortices, and baroclinity |(l/p 2 )£- (Vp x Vp)|, source 
of vorticity, at redshift z — 0. 
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Fig. 5. — The probability distribution function (PDF) of vorticity cJj at redshifts z = 4 
(square), z = 2 (triangle) and z = (cross). The solid line gives a log-normal fitting result 



We calculate the probability distribution function (PDF) of the three components of the 
vector field cJj at redshifts z — 4, 2, 0. Giving that the vorticity field is isotropic, the PDFs 
of its three components cJj, i — x,y,z should be statistically identical, which is justified in 
our samples. We take an average over these three components at these redshifts and give the 
results in Figure 5. The PDF at present epoch exhibits a long tail, and can be approximately 
fitted by a log-normal distribution as 



where the variance a = 0.98, /x = 0.37, which implies that the vorticity field is intermittent, 
i.e. the probabilities of forming big vortical structures are much larger than Gaussian fields. 

It shows that the PDF of vorticity fields has been always non-Gaussian since redshift 
z ~ 4, which is remarkably different from the velocity field of the IGM. The PDF of the 
velocity and pairwise velocity fields of dark matter and the IGM are Gaussian at high 
redshifts, corresponding to the linear phase of evolution (Yang, et al 2001). The evolution 
of the IGM vorticity field does not undergo a linear and Gaussian phase over cosmic times, 
since the vorticity can only be produced via nonlinear evolution. In this sense, the vorticity 
field is more effective than the velocity field to track the nonlinear evolution of the IGM. 

Another interesting feature indicated in Figure 5 is that the PDF at high redshifts is 
approximately of exponential, and evolves into log-normal distribution at later phase. Thus, 



for z = 0. 



4.2. PDF of the Vorticity Fields 





the PDF at different redshifts cannot be converted to each other by a scaling transformation. 
It implies that the turbulence experiences a strong nonlinear evolution, which will be revisited 
in next subsection. 



4.3. Power Spectra of Velocity and Vorticity Fields 

In a statistically homogeneous fluid, one can define the spectrum tensors $ij(k) and 
Qij (k) as the Fourier counterparts of the two-point correlation tensors of velocity (v «(x + r)vj(x)) 
and vorticity (wj(x + r)wj(x)), 

<Mk) = JL y Mx + v)v^))e-^dv (7) 

fly(k) = p y ( Wi (x + r)^(x))e- k - r rfr. (8) 

respectively, where (...) denotes average over spatial coordinates x. 
For a homogeneous turbulence, we have (Batchelor, 1959) 

%(k) = [5 l3 k 2 - kikjlQuQt) - fc 2 ^(k), (9) 

and hence, 

^(k) = fc 2 $«(k). (10) 
The power spectra of velocity and vorticity fields are defined respectively as 

p v (k) = y i^(k)5dki - k)dk-, PM = y ^ii(k)<j(iki - k)dk. (ii) 

Combining eqs. (10) and (11) yields 

P u {k) = k 2 P v (k). (12) 

This relation is an important property of homogeneous turbulence (Batchelor, 1959), and 
can be used to measure the developed level of turbulence. If the velocity and vorticity fields 
of a fluid satisfy the relation given by eq.(12), it should be in the state of fully developed 
homogeneous turbulence. Otherwise, it would be less developed. 

Figure 6 compares the power spectra P u (k) with k 2 P v (k) at z = 4 (top left), 2 (top 
right) and (bottom left), respectively. It shows that at high redshift z — 4, the power 
spectrum P u (k) is much less than k 2 P v (k) at almost all scales, which means that not all, 
actually only a small part, of the fluctuations of velocity field can be related to the random 



0.0 0.5 1.0 1.5 2.0 2.5 

lg(k/27i) (1/25h"'Mpc) 

z= 4 

R (V'Mpc) 

10.0 1.0 0.1 



16 


. i 


1 




15 
















14 








13 








12 








1 1 




lg[k 2 P„(k)] 








ig[p„M] 




10 









0.0 0.5 1.0 1.5 2.0 2.5 

lg(k/27i) (1/25h"'Mpc) 

z= 



0.0 0.5 1.0 1.5 2.0 2.5 

lg(k/27i) (1 /25h"'Mpc) 

z= 2 

R (h~'Mpc) 

10.0 1.0 0.1 



1 6 F 




1 n F , . , , , i , , , i , , = 

0.0 0.5 1.0 1.5 2.0 2.5 

lg(k/2n) (1/25h"'Mpc) 

z= 



Fig. 6. — The power spectra P UJ (k) and k 2 P v (k) at redshifts z = 4 (top left), 2 (top right ) and 
(bottom left) from 512 3 simulation. The bottom right plot gives a resolution comparison 
of these two terms at redshift z = 0. 
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field of vorticity, and the turbulence is less developed by that time. While evolving to redshift 
z ~ 2, the turbulence is developed starting from the small scale 0.2ft, -1 Mpc and up to 0.8/i" 1 
Mpc. At the present time, z — 0, the turbulence is fully developed and extended to the scale 
3/i _1 Mpc, the typical scale of a cluster. The deviations of P u (k) from k 2 P v (k) on scales 
less than 0.2 h^ 1 Mpc are probably due to the energy dissipation of turbulence to thermal 
energy, or the virialization, on small scales. A panel of these two terms in the simulation run 
of lower resolution, 256 3 , is also presented in Fig. 6 and provides a convergence test of the 
resolution effect. It shows that the resolution does affect the lower end of turbulent scale as 
a result of dissipation. However, the turbulence on large scale is resolution converged. 

Figure 6 also shows that the variance of velocity field on large scales is remarkably 
larger than that of vorticity field, especially at high redshifts. It indicates that the variance 
on large scales is not from the turbulent motion of the IGM and probably from bulk motion, 
which is due mainly to the falling into gravitational well. Therefore, to identify the variance 
of a velocity field as the signature of turbulence (e.g. Iapichino & Niemeyer 2008) may 
be questionable even on scales of clusters, as they generally contain many substructures at 
redshifts less than 2. 

We can also explore the evolution of turbulence with the spectrum of mean kinetic 
energy density E(k) defined by / °° E(k)dk = l/2(p(r)t> 2 (r)). The energy spectra E(k) at 
redshifts z = 4, 2 and are shown in Figure 7. The energy spectra can be approximately 
fitted by a power law k~ a with a — 1 in the scale range of 0.15-3 h' 1 Mpc for z = 2 and 
0.15 - 10 h^ 1 Mpc for z — 0. These scale ranges are larger than that given by the power 
spectrum of velocity and vorticity. This is probably because the turbulent flow is strong in 
high density areas. Figure 7 shows that the energy spectrum becomes very steep at scales 
less than 0.15 hr l Mpc because of the dissipation on small scales. The energy spectrum at 
z = 4 cannot be fitted with the power law of k^ 1 . It indicates that turbulence has not yet 
developed by then. Turbulence is effective at transferring kinetic energy on large scales to 
small one. Therefore, it leads to the power spectrum at z < 4 to be more flat than that of 

2 = 4. 



5. Effects of Turbulent IGM on Structure Formation 

5.1. Non-thermal Pressure 

An early attempt of including the effect of turbulent motions into gravitational collaps- 
ing processes was made by Chandrasekhar (1951). In his quantitative theory, he investigated 
the effect of micro turbulence in the subsonic regime. If turbulence is statistically homo- 
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Fig. 7. — The power spectrum of kinetic energy (solid line) at redshifts z = 4 (left), 2 (middle) 
and (right). A power law A; -1 (dashed line) is used to fit the power spectrum. 

geneous, it will contribute an extra pressure ptub = p(v 2 ) on large scales. In the linear 
regime, Chandrasekhar derived a dispersion relation by introducing an effective sound speed 
c s, e // = c s + (1/3) (v 2 ) where (v 2 ) is the rms velocity dispersion due to turbulent motion. 
Obviously, the turbulence will slow down, or even halt the gravitational collapsing. 

Chandrasekhar's result had been improved by a more elaborate investigation (Bonazzola 
et al. 1992) , in which the scale dependence of the turbulent energy was taken account in 
the analysis of system instability. Actually, the gravitational instability on a scale R will 
not be affected by fluctuation modes with wavelengths larger than R, and the fluctuation of 
velocity on the scales k < 2n/R do not contribute to the turbulent pressure for resisting on 
gravitational collapsing on scales that larger than R. Quantitatively, the turbulent pressure 
on the scale R can be estimated by (Bonazzola et al 1987) 

Ptur(k R )= / E(k)dk, (13) 
Jk R 

where E{k) is the turbulent power spectrum, = 2n/R, and k max = 2n/ldi SS is the wavenum- 
ber corresponding to the minimal scale ld iss below which the turbulence decays due to energy 
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dissipation or virialization. 

According to the results presented in §4.3, the turbulence is fully developed on scales 
from 0.2 h~ l Mpc up to 3 h~ l Mpc since redshift z ~ 2. The direct outcome of the turbulence 
on those scales is expected to alter significantly the hydrostatic equilibrium state of the IGM 
or the process of structure formation. The turbulent pressure ptur{k>R) as a function of k,R 
is shown in Figure 8, where /^ ss = 0.2h _1 Mpc inferred from the power spectrum analysis 
in §4. 3. In practical calculation, since the energy spectrum E(k) declines fast beyond 0.2 
h^ 1 Mpc (Fig. 7), one can take k max going to infinity safely. Since we have approximately 
E(k) oc k' 1 , pturiku) given by eq.(13) is weakly dependent on k. We also show the energy 
spectra E{k) in Figure 8. 

Using the power spectra measured in Figure 7, the turbulent pressure is estimated to 
be 1.5 x 10~ 17 g cm _1 s~ 2 at z — 0. According to p/p = RT/p, the effective temperature due 
to turbulent pressure is about 1.0 x 10 6 K in regions of mean overdensity and 1.0 x 10 5 K in 
regions of 10 time mean overdensity. Deduced from Lja forests of quasars, the temperature 
of IGM at p b ~ 1 — 10p bi0 is about 2 x 10 4 K. Therefore, the nonthermal pressure of the 
turbulent flow could be larger than the thermal pressure of the IGM. 




1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 

lg(k ml „) (l/25h"'Mpc) lg(><™„) (1/25h~'Mpc) 

z= 2 z= 



Fig. 8. — The spectrum of turbulent pressure Ptur(k m in), given by Ptur(k min ) = J k max E(k)dk, 
at redshifts z = 2 (left) and z = (right). The energy spectra E(k) are also shown in each 
panel. 

The scale-dependence of the turbulent pressure is very weak. A decrease of one order of 
magnitude in scales from R = 3 to 0.3 h^ 1 Mpc can only lead to deceases in the pressure pt ur 
by a factor of 4. On the other hand, the mass m of a cluster is related to its scale radius r s 
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approximately as m oc r s 3 (Cooray & Sheth, 2002). The gravitational potential of m halos at 
r s is Gm/r s oc r 2 . Therefore, the ratio of the turbulent pressure to the gravitational potential 
at r s would be larger for clusters with smaller mass m. The effect of turbulent pressure on 
gravitational collapsing of baryon gas would be more significant on smaller clusters. 




D.a 0.2 o.+ o.6 o.s i.a D.a 0.2 0.+ o.s o.s i.a 

x x 

Fig. 9. — The distribution of ln(a; 2 /2SV,SV,), which characterizes the net effect of turbulence 
on clustering and positive value represents prevention , in a 2-D slide of 25 x 25 x 0.1 h~ 3 
Mpc 3 at redshift z = 2 (left) and (right). 



5.2. Vorticity and the Growth Rate of Velocity Divergence 

In the nonlinear regime of the IGM gravitational clustering, the dynamical effect of 
turbulence can be estimated by eq.(3). Here, we are focusing on the first two terms from 
vorticity and strain rate. As discussed in §2.2, the net effect on the clustering is determined 
by the sign of quantity, 

1 1 

-u 2 - SySij = -[(diV^idiVj) - ^{djVi^diVj)]. (14) 

For a Gaussian velocity field, we have (3(djVi)(diVj)) = ((diVj)(diVj)), and in average, the 
net effect of velocity field in eq.(3) is statistically null. However, for a non-Gaussian velocity 
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field, it can be either positive or negative, which is dependent on the property of the velocity 
field. 

In homogeneous and isotropic turbulence, as ((djVi)(diVj)) = (Batchelor 1959), the 
signs of Eq.(14) are always positive, which results in prevention of gravitational collapsing 
in the IGM. Figure 9 plots the spatial distribution of ln(u 2 /2SijSij) in the same simulation 
slide as that used in Figure 3. Comparing Figure 9 with Figures 3 , we find that all those 
cells with \n(u 2 /2S ij S ij ) > are located in the clouds around density peaks, where the 
vorticity is dominant. It provides a mechanism to prevent or slow down the IGM clustering 
with respect to the underlying dark matter. 

We search for cells with (l/2)u 2 — SijSij > 0. At redshift z = 0, there is a fraction 7.6% 
of volume, 16.6% of mass, with positive values of (l/2)u 2 — S^S^, while at redshift z — 2, 
this volume fraction has decreased down to 2.6%. Thus, the effect of turbulence becomes 
stronger to prevent the IGM clustering at lower redshifts. 

z= 




□ .□1 0.10 1.00 10. DO 10D.00 1000.00 

Pc 

Fig. 10. — Comparison of effects of turbulent pressure and thermal pressure on the divergence 
d (eq.(3)), [\u 2 — SijSij]/[~V 2 p], to baryon density of randomly selected cells with \oj 2 > 
S^Sij at redshift z — 0. Solid line gives the mean value of this ratio at every density bin. 
Broken lines give the cumulative probability 20%, 50%, 70% and 90%, from bottom to up. 

In order to compare the effects of turbulent pressure and thermal pressure on the di- 
vergence d, we calculate the ratio of (l/2)w 2 — S^S^ to — V 2 p , taken from eq.(3), cell by 
cell. The result is presented in Figure 10, in which the cells are randomly selected from 
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those with ^u 2 > SijSij in the whole simulation samples. We find that for most densities 
nearly 30% of those cells with |w 2 > S^S^ have a value of this ratio larger than 1 and 
indicate that mean turbulent pressure dominates over thermal pressure in them. Obviously, 
the dynamical prevention provided by turbulence could be comparable to that of thermal 
pressure and even become dominate in a considerable fraction of the whole volume. 

6. Discussions and Concluding Remarks 

The relationship between the fields of vorticity and velocity is similar to the relationship 
between the current of charge density and magnetic field, and thus, vorticity would be a 
measurement of the coherent spatial structures of velocity field. Moreover, the dynamical 
equation of vorticity is free from the gravitational field of dark matter and cosmic fluid. 
These remarkable features are very useful to study the nonlinear behavior of cosmic baryon 
fluid, especially the clustering behavior of the turbulent cosmic fluid in the gravitational field 
of underlying dark matter. 

We show that the vorticity field of baryonic matter is significantly increasing with time 
when redshift z < 2. It can be understood that vorticity is effectively generated by shocks 
and complex structures of the baryon fluid, and then amplified by the rate-of-strain. At 
redshift z — 4, the largest vorticity is only of the order of tot ~ 10, while it is tot ~ 10 2 at 
present universe. The IGM vorticity field is non-Gaussian and intermittent at all redshifts. 
The PDF of vorticity evolves from approximately exponential distributions at high redshifts 
to a distribution with log-normal long tail at present epoch. 

The spatial configuration of the vorticity field is found to be very different from that 
of velocity and mass density. The distribution of vorticity does not follow the underlying 
matter structures, such as filaments and sheets. It always shows 3-D cloudy structures 
around gravitational collapsed regions, i.e. the knots in the filament-sheets structures. Even 
in regions surrounding high density structures, vorticity can be strong because complex 
structures, such as curved shocks and collision of shocks, are already formed around knots 
at their early phase of formation. Vorticity would be more effective to reveal the clustering 
behavior, which is overlooked by the mass density field in some way. 

The fluctuations of vorticity field is useful to measure the development of a fully devel- 
oped turbulence in the cosmic fluid. The relation between the power spectra of vorticity and 
velocity provides a measurement on the scale of velocity fluctuations where turbulence is fully 
developed. We find that the cosmic fluid is in the state of fully developed homogeneous and 
isotropic turbulence in the scale range of 0.2/i -1 Mpc to 3.0/i -1 Mpc at present epoch. With 
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this result, we calculate the turbulent pressure. It is of the order of 1.5 x 10~ 17 g cm _1 s -2 at 
z = in average, which is equivalent to the thermal pressure of gas with mean cosmic baryon 
density at temperature 1.0 x 10 6 K. It tends to slow down the gravitational clustering of the 
baryon fluid. Moreover, the spectrum of turbulence pressure is weakly dependent on scale 
k, and then the effect of turbulent pressure would be relatively stronger on smaller objects. 

The turbulent pressure may shed light on the problem of overcooling, i.e. the fraction of 
cold gas and stars in regions of galaxies, galaxies groups and clusters given by ACDM sim- 
ulations is significantly higher than the observed value at z ~ ( Nagai & Kravtsov, 2004, 
Crain et al. 2007, Keres et al. 2009). A possible way to solve this problem is to assume 
that the IGM undergo a pre-heating at low redshift (e.g. de Silva et al 2004). However, the 
pre-heating model is strongly in contradiction with the observations of the low-redshift Lja 
forest of quasars, which cannot exist if the temperature of the IGM is > 10 5 K. Galactic 
winds is another mechanism proposed to suppress star formation in galaxies. Hydrody- 
namic simulations, however, suggest that such feedback would be inefficient in galaxies with 
M ga i > 10 9 M o (Mac Low & Ferrara 1999). Turbulent pressure essentially is dynamical and 
nonthermal. It can play the similar role as thermal pressure to prevent the gravitational 
clustering, while does not affect the thermal state and ionizing process of hydrogen in the 
IGM. The turbulent IGM can remain a temperature of 10 4_5 K and hence consistent with the 
observation of Lya forest. If the IGM is turbulent, the Lya absorption lines will not only 
show thermal broadening but also turbulent broadening. Observation of Lyct line widths 
of HI and Hell indicates that the broadening of Lya forest is partially given turbulence 
broadening (Shull et al 2004, Zheng et al 2004, Liu et al 2006). 

Vorticity enhances the transportation of mass, momentum and kinetic energy. The 
cascade of vortical structures leads to transfer of kinetic energy of vortical motion from large 
scales to small scales. The turbulence energy will further dissipate into thermal motion. 
This processes will enhance efficiently the entropy production via the thermalization and 
virialization. In addition, the turbulent motion can cause diffusive mixing of materials, 
which tends to wipe out gradients in the distribution of chemical composition. The details 
will be reported in the near future. 
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A. The basic equations 

A.l. Hydrodynamical equations for the IGM 

The IGM is assumed to be an ideal fluid with polytropic index 7. The hydrodynamic 
equations for the IGM in the expanding universe can be written in the following form 

U + diMU]=f{t,U) (Al) 

where <9j = d/dXi (i = 1,2,3), denote the proper coordinates, which are related to 
comoving coordinates Xi by = a(t)xi, a(t) being the scale factor. The quantity U in 
eq.(Al) contains five components as 

U = (p,pv,E) (A2) 

where p is the comoving density of the IGM, v = {vi} (i — 1, 2, 3) are the peculiar velocity on 
three axes, E = Pj (7 — 1) + |pv 2 is the total energy per unit comoving volume, P = a 3 p, and 
p is the pressure of the IGM. The quantities fi(U) in Eq.(Al) are given by the conservation 
laws of mass, momentum and energy as 

fi(U) = [pvi,p(v!) 2 + P^pv^iPViv^v^E + P)\ 

f2(U) = [pV2,pV 1 V2,p(v2) 2 + P,pV 2 V 3 ,V 2 (E + P)] 

h{U) = [pv 3 ,pv 1 v 3 ,pv 2 v 3 ,p(v 3 ) 2 + P,v 3 (E + P)] (A3) 
The "force" term f(t, U) on the right hand side of Eq. (Al) is given by 

/(*, U) = (0, —pv + -pg, -2-E + -pw ■ g - A rad ). (A4) 
a a a a 

The term of — (a/a)pv is from the expansion of the universe. The term of A rad in Eq.(A4) is 
given by the radiative heating-cooling of the baryon gas. The gravitational force g = — V0 
is produced by the matter including CDM and baryon , given by 

V 2 = —ptoAot- (A5) 
a 

where the operator V acts on the comoving coordinate x. 5tot = [ptot( x , t) — ptot]/ptot, and 
Ptot is the total comoving mass density. Its mean value is ptotif) = oc a~ 3 . The 

gravitational potential is zero (or constant) when the density perturbation 5 tot = 0. 
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A. 2. Vorticity equation 

Euler equation in comoving coordinates 

d t pvi + -(djpvjVi + djSijP) = --pvi + -pgi, (A6) 

(X (X CL 

or 1 '1 

pd t Vi + Vid t p + -(vidjpvi + pv ] djVi + d^P) = — pv { + -pg h (A7) 

(X (X (X 

, where d t = d/dt. Therefore 

1 (X 1 

pd t Vi + -(pVjdjVi + djSijP) = — pvi + -pg h (A8) 

(X (X (X 

111 

dj5ij Pd t Vi + -VjdjVi = — {-d^yP + - (A9) 
a a p 

Using Levi Civita symbol 

tijkUlm = 8jl8km — Sj m Skl (AlO) 

we have 

V jdj V i = -jdi V j V j - e ijkVjOJ k , (All) 

where u>i = €ij k d^v k is vorticity. 

Taking operator of curl e^ k & on eq.(A9), we have term by term. 

e ijk d J d t v k = d t e ijk &>v k = d t Ui (A12) 
tijkd 3 vidiv k = e ijk d j ^d k vivi - e ijk d J e klm viOj m (A13) 

£ijkd J e k i m viu m = e kij e klm d j viuj m = d m ViU rn - d l viUi = u m d m Vi - u^vi - vid l Ui (A15) 



e ijk d j -d k v lVl = (A14) 



Therefore, we have vorticity equation as 

d t Ui + -vid l Ui = -(u m d m Vi - Uid l vi + \e ijk djpd k p - awi). (A16) 
a a p z 

Because Ui(d j Vi — d l Vj) = 0, we have 

11 1 

d t Ui + ~vid l Ui = -(SijUj - dui + -^e ijk djpd k p - auji) (A17) 
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where 

S ij = -(&v i + &v j ) (A18) 

and d = diVi. In vector format 

11 1 
d t w + -v-Vu = -(S • w - rfw + -^Vp xVp - aw) (A19) 
a a p^ 



A. 3. Equation of divergence 

Taking operator d l on eq.(A9), we have term by term, 

tfdtVi = d t d (A20) 
d l Vjd,jVi = Vjdjd l Vi + (d t v j )(d j v i ) = vjdjd + (d t v j )(d j v i ) (A21) 

Using 

& Vj = l - (d\ + d j v i ) + ^ (&Vj - d j Vi ), ( A22) 

we have 

1 ...... . 1 

{& , Vj){djV i ) = SijSij + -{d'vj - WvijicPvi - d l Vj) = S^S^ + -e ijk djV k ea m div m . (A23) 

Therefore, the equation of divergence is 



d t d + -Vi&d = -(^UiUi - SyS^ - -didip + ^djpdjp -ad - -^-(p - p ))- (A24) 
or in vector format 

d t d + -v • Vd = -{\u ■ to - SijS tJ - -V 2 p + ^(Vp) • (Vp) -ad- — (p - p )). (A25) 
a a 2 p p z a 



A. 4. A brief description of the numerical algorithm. 

We use the fifth order finite difference WENO scheme (Jiang & Shu 2006) to demonstrate 
the basic idea of the WENO methodology. The fifth order WENO finite difference spatial 
discretization to a conservation law such as 

u t + f(u) x + g{u) y + h{u) z = (A26) 
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approximates the derivatives, for example f(u) x , by a conservative difference 

f(u) x \ x=Xj « ^- (/ j+1/2 - /,_ 1/2 ) (A27) 

along the x axis, with y and z fixed, where /j+1/2 is the numerical flux. <?('u)y and h{u) z 
are approximated in the same way. Hence finite difference methods have the same format 
for one and several space dimensions, which is a major advantage. For the simplest case of 
a scalar equation (A26) and if f'(u) > 0, the fifth order finite difference WENO scheme has 
the flux given by 

A+i/2 = wJ^ 1/2 + wJ^ 1/2 + wJf} l/2 (A28) 
where are three third order accurate fluxes on three different stencils given by 

= \flVi-*) - \fiuj-,) + y /(«,-), (A29) 
/S/2 = 4/K-i) + f /(«;) + (A30) 



/S/2 = + - ^/(«i + 2). (A31) 



Notice that the combined stencil for the flux /j+1/2 is biased to the left, which is upwinding 
for the positive wind direction due to the assumption f'(u) > 0. The key ingredient for the 
success of WENO scheme relies on the design of the nonlinear weights w i: which are given 
by 

Wi = ^3 w k = 7— 0^2, (A32) 

where the linear weights 7^ are chosen to yield fifth order accuracy when combining three 
third order accurate fluxes, and are given by 

7l = To' 72 = 5' 73 = Id ; (A33) 



the smoothness indicators (3k are given by 

Pi = ^ (/K-2) - 2/( Uj _i) + fiuj)) 2 + ± (/(uj-a) - 4/( Uj -_i) + 3/(^)) 2 (A34) 

13 1 

& = ^ (/(uj-i) - 2f( Uj ) + f{u j+1 )f + i (/(«,-_!) - /(^ +1 )) 2 (A35) 

A = ^ (/(%) - 2/(u j+ i) + /(% +2 )) 2 + ^ (3/(uj) - 4/(u j+1 ) + /(u j+2 )) 2 , (A36) 

and they measure how smooth the approximation based on a specific stencil is in the target 
cell. Finally, e is a parameter to avoid the denominator to become and is usually taken 
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as e = 10 6 in the computation. There are no other parameters needed to be tuned by the 
user in the WENO method. 

Meanwhile, the time step in simulation is set to the minimum value among two time 
scales. One is given by Courant condition as 

< cflwm (A37) 

max{\vi + c s , v 2 + c s , v 3 + c s ) 

where Ax is the cell size, c s is the local sound speed, vi, v 2 , and v 3 are fluid velocities, and 
CFL is the Courant number, here CFL = 0.60. The other one is from cosmic expansion, 
which requires that Aa/a < 0.02 within a single time step. 

The WENO scheme is proven to be uniformly fifth order accurate including at smooth 
extrema, and this is verified numerically. Near discontinuities the scheme produces sharp 
and non-oscillatory discontinuity transition. The approximation is self-similar. Namely, 
when fully discretized with the Runge-Kutta methods, the scheme is invariant when the 
spatial and time variables are scaled by the same factor. This is a major advantage for 
approximating conservation laws which are invariant under such scaling. 
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